// $Id: sPHENIX_MCApplication.cxx 536 2011-03-09 17:06:04Z ivana $

#include "sPHENIX_MCApplication.h"
#include "sPHENIX_MCStack.h"
#include "sPHENIX_MagField.h"

#include <TROOT.h>
#include <TInterpreter.h>
#include <TVirtualMC.h>
#include <TPDGCode.h>
#include <TGeoManager.h>
#include <TVirtualGeoTrack.h>
#include <Riostream.h> 
#include <TCanvas.h>
#include <TParticle.h>

#include <sys/time.h>
#include <stdio.h>
#include <unistd.h>

/// \cond CLASSIMP
ClassImp(sPHENIX_MCApplication)
/// \endcond

//_____________________________________________________________________________
sPHENIX_MCApplication::sPHENIX_MCApplication(const char *name, const char *title, Int_t nOfEvents, FileMode fileMode)
: TVirtualMCApplication(name,title),
  Stack(0),
  DetConstruction(0),
  EMCalSD(0),
  GEMSD(0),
  PreshowerSD(0),
  HCalSD(0),
  RootManager("sphenix", fileMode),
  MagField(),
  NofEvents(nOfEvents)
{
  /// Standard constructor
  /// \param name   The MC application name 
  /// \param title  The MC application description
  /// \param fileMode  Option for opening Root file (read or write mode)
  
  // Create a user stack
  Stack = new sPHENIX_MCStack(100); 
  
  // Create detector construction
  DetConstruction = new sPHENIX_DetectorConstruction();
  
  // Create EMCal SD
  EMCalSD = new sPHENIX_EMCalSD("EMCal",DetConstruction);

  // Create a GEM SD
  GEMSD = new sPHENIX_GEMSD("GEM",DetConstruction);
  
  // Create a Preshower SD
  PreshowerSD = new sPHENIX_PreshowerSD("Preshower",DetConstruction);
  
  // Create a HCal SD
  HCalSD = new sPHENIX_HCalSD("HCal",DetConstruction);

  MagField = new sPHENIX_MagField();

  // Give RootManager pointers
  
  RootManager.SetDetectorConstruction(DetConstruction);
  RootManager.SetEMCalSD(EMCalSD);
  RootManager.SetPreshowerSD(PreshowerSD);
  RootManager.SetHCalSD(HCalSD);
  RootManager.SetNofEvents(NofEvents);
  RootManager.InitializeHistograms();

  // Constant magnetic field (in kiloGauss)
  //MagField = new sPHENIX_MagField(20., 0., 0.);
  
  // It si also possible to use TGeoUniformMagField class:
  // MagField = new TGeoUniformMagField(20., 0., 0.);
}

//_____________________________________________________________________________
sPHENIX_MCApplication::sPHENIX_MCApplication()
  : TVirtualMCApplication(),
    Stack(0),
    DetConstruction(0),
    EMCalSD(0),
    GEMSD(0),
    PreshowerSD(0),
    HCalSD(0),
    NofEvents(0)
{    
/// Default constructor
}

//_____________________________________________________________________________
sPHENIX_MCApplication::~sPHENIX_MCApplication() 
{
  /// Destructor  
  
  delete Stack;
  delete DetConstruction;
  delete EMCalSD;
  delete GEMSD;
  delete PreshowerSD;
  delete HCalSD;
  delete gMC;
  gMC = 0;
}

//
// private methods
//
//_____________________________________________________________________________
void sPHENIX_MCApplication::SetInputFile(const char *in_file) {
  input_file = string(in_file);
  use_input_file=1;
  cout << "Opening up file " << input_file.c_str() << endl;
  TFile *f1 = new TFile(input_file.c_str());
  input_tree = (TTree*)f1->Get("t1");
  input_tree->SetBranchAddress("particlelist",&input_part_list_ptr);
  event_number = 0;
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::RegisterStack()
{
  // Register stack in the Root manager.
    
  RootManager.Register("stack", "sPHENIX_MCStack", &Stack);   
}  

//
// public methods
//

//_____________________________________________________________________________
void sPHENIX_MCApplication::InitMC(const char* setup)
{    
/// Initialize MC.
/// The selection of the concrete MC is done in the macro.
/// \param setup The name of the configuration macro 

  gROOT->LoadMacro(setup);
  gInterpreter->ProcessLine("Config()");
  
  gMC->SetStack(Stack);
  gMC->SetMagField(MagField);
  gMC->Init();
  gMC->BuildPhysics(); 
  
  //RegisterStack();
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::RunMC()
{    
/// Run MC.
/// \param nofEvents Number of events to be processed

  gMC->ProcessRun(NofEvents);
  FinishRun();
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::FinishRun()
{    
/// Finish MC run.
  TCanvas *can1 = new TCanvas("can1","can1",0,0,1000,1000);
  gGeoManager->GetMasterVolume()->Draw();
  gGeoManager->SetVisOption(0);
  gGeoManager->SetTopVisible();
  gGeoManager->SetVisLevel(0); 
  gGeoManager->DrawTracks("/*");

  RootManager.WriteAll();
  cout << " DONE " << endl;
  gROOT->GetListOfCanvases()->ls();
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::ConstructGeometry()
{    
/// Construct geometry using detector contruction class.
/// The detector contruction class is using TGeo functions or
/// TVirtualMC functions (if oldGeometry is selected)

  // Cannot use Root geometry if not supported with 
  // selected MC
  cout << "Geometry will be defined via TGeo" << endl;
  DetConstruction->ConstructMaterials();  
  DetConstruction->ConstructGeometry(); 
} 

//_____________________________________________________________________________
void sPHENIX_MCApplication::InitGeometry()
{    
/// Initialize geometry
  
  EMCalSD->Initialize();
  GEMSD->Initialize();
  PreshowerSD->Initialize();
  HCalSD->Initialize();
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::GeneratePrimaries()
{    
/// Fill the user stack (derived from TVirtualMCStack) with primary particles.
  
 // Track ID (filled by stack)
 Int_t ntr;
 
 // Option: to be tracked
 Int_t toBeDone = 1; 
 
 // Particle type
 //Int_t pdg  = kElectron;
 Int_t pdg = 211;

 // Polarization
 Double_t polx = 0.; 
 Double_t poly = 0.; 
 Double_t polz = 0.; 

 // Position
 Double_t vx  = 80.; 
 Double_t vy  = 80.; 
 Double_t vz =  -50.;
 //Double_t vz = -0.5*(DetConstruction.GetWorldFullLength());
 Double_t tof = 0.;

 // Energy
 Double_t kinEnergy = 10.;
 // Double_t mass = 0.000511;  // elek
 Double_t mass = 0.135;
 Double_t e  = mass + kinEnergy;
 
 // Momentum
 Double_t px, py, pz;
 px = 0; 
 py = 0; 
 pz = sqrt(e*e - mass*mass); 

  // Add particle to stack 
  if (!use_input_file) {
   Stack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, vx, vy, vz, tof, polx, poly, polz, kPPrimary, ntr, 1., 0);
  } else {
   input_tree->GetEntry(event_number);   for (int i=0; i<input_part_list_ptr->size(); i++) {
     TParticle *temp = input_part_list_ptr->at(i);
     //cout << "Particle " << i << " has energy = " << temp->Energy() << endl;
     if (temp->GetStatusCode() == 1) {
       pdg = temp->GetPdgCode();
       px = temp->Px(); py = temp->Py(); pz = temp->Pz(); e = temp->Energy();
       vx = temp->Vx()/10.; vy = temp->Vy()/10.; vz = temp->Vz()/10.; tof = temp->T();
       Stack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, vx, vy, vz, tof, polx, poly, polz, kPPrimary, ntr, 1., 0);
       ntr++;
     }
   }
   event_number++;
 }

}

//_____________________________________________________________________________
void sPHENIX_MCApplication::BeginEvent()
{    
/// User actions at beginning of event.
/// Nothing to be done this example

  Stack->SetObjectNumber();
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::BeginPrimary()
{    
/// User actions at beginning of a primary track.
/// Nothing to be done this example
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::PreTrack()
{    
/// User actions at beginning of each track.
/// Nothing to be done this example
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::Stepping()
{    
/// User actions at each step

  EMCalSD->ProcessHits();
  GEMSD->ProcessHits();
  PreshowerSD->ProcessHits();
  HCalSD->ProcessHits();
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::PostTrack()
{    
/// User actions at each step.
/// Nothing to be done this example
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::FinishPrimary()
{    
/// User actions after finishing of a primary track.
/// Nothing to be done this example
}

//_____________________________________________________________________________
void sPHENIX_MCApplication::FinishEvent()
{    
/// User actions after finishing of an event
/// Nothing to be done this example

  // Geant4 own visualization is activated via G4 macro (g4config.in)
 
  // TGeo visualization
  if ( gGeoManager && 
       gGeoManager->GetListOfTracks() &&
       gGeoManager->GetTrack(0) &&
       ((TVirtualGeoTrack*)gGeoManager->GetTrack(0))->HasPoints() ) {
       
     gGeoManager->SetVisOption(0);     
     gGeoManager->SetTopVisible();
     gGeoManager->DrawTracks("/*");  // this means all tracks
  }    

  //  EMCalSD->Print();
  //  GEMSD->Print();

  RootManager.Fill();
  RootManager.RunAnalysis();

  EMCalSD->EndOfEvent();
  GEMSD->EndOfEvent();
  PreshowerSD->EndOfEvent();
  HCalSD->EndOfEvent();

  //Stack->Print();  
  Stack->Reset();
} 

//_____________________________________________________________________________
void  sPHENIX_MCApplication::ReadEvent(Int_t i) 
{
/// Read \em i -th event and prints hits.
/// \param i The number of event to be read    
  
  EMCalSD->Register();
  GEMSD->Register();
  PreshowerSD->Register();
  HCalSD->Register();

  RegisterStack();
  RootManager.ReadEvent(i);

  //Stack->Print();  
  //EMCalSD->Print();
  //GEMSD->Print();
}
